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ABSTRACT 

We perform 3-D dust radiative transfer (RT) calculations on hydrodynamic simulations of isolated 
and merging disk galaxies in order to quantitatively study the dependence of observed-frame sub- 
millimeter (submm) flux density on galaxy properties. We find that submm flux density and star 
formation rate (SFR) are related in dramatically different ways for quiescently star-forming galax- 
ies and starbursts. Because the stars formed in the merger-induced starburst do not dominate the 
bolometric luminosity and the rapid drop in dust mass and more compact geometry cause a sharp 
increase in dust temperature during the burst, starbursts are very inefficient at boosting submm flux 
density (e.g., a > 16x boost in SFR yields a < 2x boost in submm flux density). Moreover, the ratio 
of submm flux density to SFR differs significantly between the two modes; thus one cannot assume 
that the galaxies with highest submm flux density are necessarily those with the highest bolometric 
luminosity or SFR. These results have important consequences for the bright submillimeter-selected 
galaxy (SMG) population. Among them are: 1. The SMC population is heterogeneous. In addition 
to merger-driven starbursts, there is a sub-population of galaxy pairs, where two disks undergoing a 
major merger but not yet strongly interacting are blended into one submm source because of the 
large (> 15", or ~ 130 kpc at z = 2) beam of single-dish submm telescopes. 2. SMGs must be very 
massive (M* >6x 10 10 M Q ). 3. The infall phase makes the SMG duty cycle a factor of a few greater 
than what is expected for a merger-driven starburst. Finally, we provide fitting functions for SCUBA 
and AzTEC submm flux densities as a function of SFR and dust mass and bolometric luminosity 
and dust mass; these should be useful for calculating submm flux density in semi-analytic models and 
cosmological simulations when performing full RT is computationally not feasible. 
Subject headings: galaxies: high-redshift — galaxies: interactions — galaxies: starburst — infrared: 
galaxies — radiative transfer — submillimeter: galaxies 



1. INTRODUCTION 



Su bmillimeter-selected galaxies (SMGs; iSmail et"al 



19971 : iBarger et alj|1998t iHughes et alJll998tlEales et al 



1999 : see iBlain et all 120021 for a review) are extremely 



lumin ous (bolometric lu minosity Lboi ~ 10 — 10 Lp>\ 
e.g.. iKovacs et aT]|2006() . high-redshift ()Chapman et al.l 
120051 ) galaxies power ed primarily by star formation 



rather than AGN 
Valiante et all J2007t 



I Alexander et al 



y star tor 
1 l2005airbl 



2008; rbr 



2007 



_ Mcnendcz-Dclmcstre et al _ 
20091: IPope et al.1 120081: lYounger et al.l 120081 I2009bl) 
Because of their high dust content, SMGs emit almost 
all of their luminosity in the IR. As the name suggests, 
a galaxy is defined as an SMG if it is detected in the 
submm (historically, 850 /xm flux density S^o > 3 — 5 
mJy; the nature of the population is sensitive to the 
adopted flux density cut, so we define an SMG as a 
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source with 5ssn > 3 mJy), which require s Lir > 10 12 L Q 
(jKovacs et alJ l2006t iCoppin et all 120081 ) . so SMGs are 
typically ultra-luminous infrared galaxies (ULIRGs). 
Locally, ULIRGs are almos t exclusively merg i ng ga laxies 
(jSanders fc Mirabell 119961 : lLonsdale et~aTT 120061) . so 
one might expect that at least some SMGs are also 
merging galaxies. Indeed, many o bservations s upport 
a mer ger o ri gin for SMGs (e.g., Ilyison et al 
2007 1 120101 : IChapman et al.l 120031 : iNeri et al ' 
Smail et al.l 1200^ 1 ' 
20051: 



Swinbank et all 12004 
Tacconi et alJ 120061 120081: iBouche et al 





2002 




2003 



Greve et all 



Biggs fc Ivisonl |20Q. 



_ 2007; 
Capak et all 120081: lYounger et al.1 



2008L 120101: llono et all 120091: lEngel et all 120101) . Fur- 
thermore, in iNaravanan et all (|2010aL hereafter N10) 
we combined hydrodynamic simulations and radia- 
tive transfer (RT) calculations to show that major 
mergers can reproduce the full range of submm flux 
densities and typical UV-mm spectral energy distri- 
bution (SEP) of SMGs (cf. IChakr abarti et al.l [2001 
IChakrabarti fc Whitney! I2009P ). Semi-analytic models 
also predict that the SMG population is dominated 
by merger-induce d starbursts rather than quiescent 
star formation dBaugh et al.l 120051: iFontanot et all 



20071: iSwinbank et al.1 120081: lLo Faro et all 120091: 
Fontanot fc Monacol 1201(1 iGonzalez et alJl201li but cf. 



branato et"aU f2004 

However, because of the much greater rate of gas 
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supp l y onto galaxies at high redshift ( e.g.. iKeres et all 
20051 lDekel"etaLll2009l) gas fractions (lErb et all 120061: 



Tacconi et al.ll2006l 120 1 0: ibaddi et al.ll2010D andstarfo r- 
mation rates (paddi et all 120071 : iNoeske et all l2007al |O) 
of galaxies at fixed galaxy mass increase rapidly with 
redshift. Thus, at z ~ 2 — 3 even a "normal" star- 
forming galaxy can r each ULIRG luminositie s (e.g. , 
Hopkins et all l2008cl \20Wt . IDaddi et all [2001 120071: 



Dannerbauer et al. 1 120091) . Furthermore, roughly esti- 



mating submm counts using estimates of high-redshift 
major merger rates and the short duty cycle of merger- 
induced starbursts suggests that there may not be 
enoug h major mergers to account for the SMG popu- 
lation (|Dave et al.ll2010l ). This motivates the view that, 
instead, typical SMGs may be massive, gas-rich disks 
quiescently forming stars and fueled by co ntinuous gas 
supp ly from mergers and smooth accretion (jCarilli et all 
[2010. but cf. IDaddi et al.ll2009bD . 

The mode of star formation responsible for the major- 
ity of the SMG population is still a matter of debate, as 
it is difficult to discriminate between the two scenarios 
given the currently available data. A better understand- 
ing of the submm galaxy selection can clarify the nature 
of the SMG population. 

Since SMGs have redshifts z ~ 1 — 4 



(iDannerbauer et all l2002t iChapman et al. 




2005; 


Younger et all 20071 120081: Caoak et al. 


2008; 


Greve et al. 20081; ISchinnerer et all 120081; IDaddi et al. 


2009allbl; IKnudsen et al.H2010D, the observed submm flux 



density traces rest-frame ~ 150 — 400 fim, longward of 
the peak of the IR SED. Thus the observed submm flux 
density is sensitive to both the total IR luminosity and 
the "dust temperature"0 of the SED, which depend on 
the luminosity from stars and AGN absorbed by the 
dust, the mass and composition of the dust, and the 
spatial distribution of stars, AGN, and dust. Galaxies 
do not have identical SED shapes, so the dependence on 
dust temperature implies that galaxies with the highest 
submm flux density are not necessarily those with the 
highest bolometric luminosity. Furthermore, because 
star formation histories are more complicated than an 
instantaneous burst, the luminosity and instantaneous 
SFR are not necessarily linearly proportional. Thus 
the relationship between submm flux density and SFR 
is potentially more complicated than the relationship 
between submm flux density and bolometric luminosity. 
We therefore cannot say a priori that the galaxies with 
the highest submm flux densities are the most rapidly 
star-forming or most luminous bolometrically. Indeed, 
it has already been observationally demonstrated that 
submm selection does not select all the brightest galaxies 
in a given volume, as there are galaxies with luminosities 
and redshifts comparable to those of SMGs that are 
undetected in the submm because of their relatively 



hot SEDs (IChapman et al.l 120041: IChapman et al, 
Casey et al.l 120091;" Casey et al.l 120101 ; iHwang et al 



2010; 



2010; 



Magdi7eOJl2010; M agnelli et alT l2010f). A submm 



galaxy selection is clearly biased toward cold galaxies; 

8 As is convention, we will use the term "dust temperature" to 
denote the temperature derived from a single- temperature modified 
blackbody fit to the SED. This is simply a parameterization of the 
SED shape rather than a physical temperature. In our simulations 
dust grains have a continuum of temperatures, depending on both 
grain size and the local radiation field heating the dust. 



however, the details of the selection bias are yet to be 
understood. 

Despite the basic physical reasons that one does not 
expect a simple relation between submm flux density 
and SFR, a linear relation between submm flux and 
SFR has been used explicitly (and, even more fre- 
quently, implicitly) to infer SFR from observed submm 
flux densities (e.g . . IChapman et al.l 120001; iPeacock et aT 



2000; iBlain et all 120021; iScott et all l2002t iWebb et al 
20031: Ivan Kampen et al.l 120051 iTacconi et al.l 120081 : 
iWang et all 120 111 ) , typically because the data sets do 
not have enough photometric data points to precisely 
constrain the IR SED shape (Her schel data are already 
helping greatly in this regard; e.g., IChapman et al.|[2010t 



IDannerbauer et all l2010t iMagnelli et al.ll2010D. Fu rther- 
more, some theoretical studies (jDave et all 120101 have 
assumed that SMGs are the most rapidly star-forming 
galaxies in order to identify SMGs in cosmological simu- 
lations without performing RT. If SFR and submm flux 
density are not simply related this approach is problem- 
atic. 

It is clear that a better understanding of the relation- 
ship between submm flux density and SFR and, more 
generally, what galaxy properties a submm galaxy se- 
lection selects for, is needed. In other work we have 
combined hydrodynamic simulations and dust RT to 
show that major mergers of massive, gas-rich disk galax- 
ies can repro duce the 850 um flux densities (N10), 
CO properties (|Narayanan et al.ll2009f ). number densities 
(jHavward et alll2011l C. Havward et al. 2011, in prepa- 
ration), and intersec tion with the dust-obscu red galaxy 
(DOG) population (jNaravanan et al.l l2010bl) of SMGs. 
Motivated by the success of our simulations in repro- 
ducing a variety of SMG properties, we use them here 
to quantify how submm flux density depends on SFR, 
Lboi, dust content, and geometry. The aim of this study 
is to clarify for what galaxy properties a submm selec- 
tion criterion selects and to provide a discriminant among 
the different modes of star formation that could power 
SMGs. 



We 



2. METHODS 
combine hig h-resolut ion Gadget- 



2 (|Springel et al.l 12001 ISprineell 120051) 3-D N- 
body/smoothed-particle h ydrodynamics (SPH) sim- 
ulatio ns with the Sunrise (|Jonssonll2006l : Uonsson et all 
l2010f ) polychromatic Monte Carlo dust RT code in 
order to predict the submillimeter flux densities of 
high-redshift isolated and merging disk galaxies. The 
simulations presented here, part of a larger suite to 
be presented in C. Hay ward et al. (2011, in prepa- 
ration), are described in N10, so here we will only 
summarize and describe differences from that work. 
This combination of Gadget-2 and Sunrise has been 



SINGS dKennicutt et al.H2003HDale et all 


2007) galaxies 


(Jonsson et all 20101); local ULIRGs ( 


Younger et al. 



2009a ); mass i ve, q ui escen t, compact z ~ 2 galaxies 
(iWuvts et all 120091 I2010D : 24 ^m-selected galaxies 



^Narayanan et al. || 2010bf) : K+A / post -starburst galaxie s 



(jSnvder et al.ll2011f) : and XUV disks (|Bush et alll20Tol) . 
among other populations. The success of our approach 
at modeling diverse galaxy populations — both local 
and high-redshift — lends credibility to its application to 
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modeling SMGs. 

2.1. Hydrodynamic simulations 

Gadget- ^1 (ISpringel et all f200H ISpringel l200l is 
a TreeSPH (jHernquist fc Katzl 11983) code that com- 
putes gr avitational interacti ons via a hierarchical tree 
method (IBarnes fc Hut]|1986D and gas dyn amics via SPH 
(jLucvi [l977t iGingold fc Monaghanlll977l). It conserve s 
both energy and entropy (|Springel &: Hernquistl 2002) . 
The s imulations i nclude radiative heating and cooling 
as in iKatz et all (|1996l) . Star formation is modeled 
via the volume-densit y dependent Kennicutt-Schmidt 
law (|K cnmcu tt 1998a), Psfr oc pl±, with a minimum 
density threshold; this i ndex is consistent with obser- 
vations of z ~ 2 disks dKrumholz fc Thompson! I2007t 
iNaravanan et all l2008al |2011|). The density threshold 
used is n ~ 0.1 cm -3 , much less than that of the 
dense molecular gas from which stars form (n ~ 10 2 — 
10 3 cm -3 ). For this reason, and because we do not 
track the formation of molecular gas, the KS law em- 
ployed should be considered simply an empirically- and 
physically-motivated prescription to encapsulate physics 
we do not resolve. The SF prescription has been cal- 
ibrated to reproduce the global K-S law (see §3 of 
ISpringel &: Hernquistl 120 03;). Recently, some authors 



(e.g., iHopkins et al.l 120111 ) have presented simulations 



that resolve the density threshold for molecular gas for 
mation; we plan to compare such simulations to our cur- 
rent simulations in future work. 

The structure of the ISM is modeled via a two-phase 
sub-resolution model in which cold , dense clouds are em- 
bedd ed in a diffuse, hot medium ([Springel &: Hernquistl 
2003). This medium is pressurized by supernova feed- 
back, which heats the diffuse ISM and evaporates the 
cold clouds ()Cox et al.ll2006bT ). Metal enrichment is cal- 
culated by assuming each gas particle behaves as a closed 
box. Black hole particles accrete via Eddington-limited 
Bondi-Hoyle accretion and deposit 5% of their emitted 
luminosity — calculated from the accretion rate assum- 
ing 10% radiative efficiency, Lboi = O.lMc 2 — to the sur- 
rounding ISM as thermal energy (ISpringel et al.1 120051 : 
iMatteo et al.ll2005l ) . We refer the reader to the references 
given above for the full details of the Gadget-2 code and 
the sub-resolution models employed. 

We focus on two simulations, one isolated disk and one 
major merger. We embed exponential disks with bary- 
onic mass 4 x 10 n Mp, in 9 x 10 12 M^ dark matter ha- 
los described by a iHernquistl (|1990f l density profile. The 
disks are i nitially 60% gas and are scaled to z ~ 3 as de- 
scribed in iRobertson et all (2006a. b). The gravitational 
softening lengths are 200/i -1 pc for the dark matter par- 
ticles and 100/i -1 pc for the star, gas, and black hole par- 
ticles. We use 6 x 10 4 dark matter, 4 x 10 4 stellar, 4 x 10 4 
gas, and 1 black hole particle per disk galaxy. For the ma- 
jor merger, we initialize two such disks on parabolic or- 
bits with initial separation i?i n it = 5i? v j r /8 and pericen- 



tric dis tance twice the disk scale length (jRobertson et al 
2006b). The orbit we focus on is the 'e' orbit o fC ox et al 
(2006a). We have checked that the differences between 
quiescent star formation and starbursts are insensitive 



to orbit as long as a strong starburst is induced (some 
orbits do not induce strong starbursts, but those are irrel- 
evant for the purpose of studying the differences between 
starbursts and quiescent star formation), and the larger 
suite of simulations used to derive the fitting functions 
includes a variety of orbits. 

2.2. Radiative transfer 

Every 10 Myr we save snapshots of the Gadget- 
2 simulations an d use the 3-D Monte Carlo dust RT 
code SunriseH (| Jonssonl [20061 : Uonsson et al.l l20ldh in 
post-processing to calculate the SEDs of the simulated 
galaxies. While we will summarize the key fe atures of 
Sunr ise h ere, we encourage th e reader to see Uonssonl 
(2006) and Uonsson et all (|2010D for full details. Except 
where noted, we us e the fiducial parameters given in 
Uonsson et al.1 (|2010f ). Sunrise calculates the emission 
from the stars and AGN in the Gadget-2 simulations 
and the a ttenuation and re-em ission from dust. Star- 
burst99 (jLeitherer et al.|[l999t) SEDs are assigned to all 
star particles according to their ages and metallicitics. 
Star particles present at the start of the Gadget-2 sim- 
ulation are assigned ages assuming that their stellar mass 
was formed at a constant rate equal to the star formation 
rate of the initial snapshot and gas and stellar metallic- 
ities via a closed-box model, Z = —yhxf g , where f g is 
the initial gas fraction and y — 0.02. Black hole parti- 
cles are assign e d the luminosity-dependent templates of 
IHopkins et all (pool by assuming that the bolometric 
luminosity of a black hole particle is Lboi = O.lAfc 2 , 
where M is the black hole accretion rate from the 
Gadget-2 simulations. 

To calculate the dust density, and thus optical 
depth along a given line-of-sight, Sunrise projects the 
Gadget-2 gas-phase metal density onto a 3-D adaptive 
mesh refinement grid using the SPH smoothing kernel. 
We have assumed 40% of the metals are in dust (jDwekl 
119981 : Uames et alJl2002j ) . We use a maximum refinement 
level of 10, resulting in a minimum cell size of 55/i _1 
pc. By performing runs with higher levels of refinement 
we have checked that the observed-frame submm flux 
density is converged to within 10%. We use the Milky 
Way R=3.1 du st model of iWeingartner fc Draini (|2001[1 
as updated by iDraine &: Lil (|2007l) . Dust models with 
different FIR opacity will lead to different relationships 
between submm flux density and dust mass, but we show 
how to r escale for different values of the opacity in Equa- 
tion |A9] 

Once the dust density grid is constructed and the in- 
put sources are assigned SEDs, Sunrise performs Monte 
Carlo RT by randomly emitting photon packets from the 
source particles and randomly drawing interaction opti- 
cal depths from the appropriate probability distribution. 
We use 10 7 photon packets total for each stage of the 
RT, having confirmed that this results in Monte Carlo 
noise of less than a few percent. The photon packets are 
scattered and absorbed by dust as they traverse the ISM. 
For each grid cell, the temperature of each dust species is 
calculated assuming the dust is in thermal equilibrium, 
and the dust re-emits the absorbed energy as a modified 
blackbody. 



9 A public version of Gadget-2 is available at 10 Sunrise is publicly 

http: / /www.mpa-garching .mpg . de/$\sim$volker/gadget/index .html http : //code .google . com/p/sunrise/ 



available 



at 
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A Sunrise feature key to this work is its treatment 
of dust self-absorption. In high-density regions, the dust 
can be opaque to its own emission, so the contribution 
of the dust emission to dust heating must be computed 
in addition to the contribution from stars and AGN. 
Sunrise computes the equilibrium dust temperatures 
self-consistently by iteratively performing the transfer of 
the dust emission and the temperature calculat ion us- 
ing a reference field technique similar to that of Uuvelal 
(2005) . (The details of the Sunrise implementatio n 
are in lJonsson et al.ll2010l and lJonsson fc Primack 2010.) 
This algorithm ensures accurate dust temperatures, and 
thus submm flux densities, even for the highly optically 
thick central starbursts. 

The results of the Sunrise calculation are spatially 
resolved, multi- wavelength (we use only 120 wavelengths 
here because of memory constraints) SEDs observed from 
7 cameras distributed isotropically in solid angle, though 
in this paper we onl y utilize the integrat ed flux densi- 
ties in the SCUBA (jHolland et al.l 11999ft and AzTEC 
(jWilson et al.l 120081 ) bands. For the purpose of calcu- 
lating observed flux densities we assume the simulated 
galaxies are at redshift z = 2. 

2.2.1. Differences from Narayanan et al. 

The primary d ifference between our sim ulations and 
those of N10 and iNaravanan et al.l (|2010bl ) is the treat- 
ment of the ISM on sub-resolution scales. In order to 
model the effects of HII and photodissociation regions 
(PDRs), Sunrise assigns star particles with ages less 
th an 10 Myr SEDs fro m the HII region template library 
of iGroves etU] (f2008ft . The time -averaged fraction of 
solid angle obscured by the PDR, /pdr, strongly affects 
the resulting a ttenuation and IR emission (for a detailed 
discussion see [Groves et all 120081 ). Narayanan et al. as- 
sumed /pdr = 1 (so that the young stars are completely 
obscured by PDRs for 10 Myr) in order to match the 
observed range of 850 /xm flux densities. Furthermore, 
they neglected the dust associated with the cold phase 
of the lSpringel fc Hernquisti (|2003| ) ISM model, typically 
> 30% of the total gas mass and > 90% of the gas mass 
in the central regions for snapshots classified as SMGs. 

Motivated by c oncerns over applicability of the 
IGroves et al.l (|2008|) models to the extreme ISM densi- 
ties and pressures encountered in our simulations, we set 
/pdr = 0, eliminating all significant dust obscuration 
from the sub-resolution PDR model. Instead, we use the 
total gas density in the SPH simulations (i.e., both the 
diffuse and cold phases) to calcula te the dust density . 
Since the dust mass implicit in the IGroves et al.l (2008) 
PDRs is not tied to the total dust mass of the simulated 
galaxy, it is possible that one can have more dust mass in 
the sub-resolution PDRs than is available in the galaxy. 
It is also possible that the sum of the dust mass in the 
PDRs is less than the total available in the cold phase 
of the sub-resolution ISM. Our treatment ensures that 
neither scenario can occur. 

Our assumed ISM structure (cold phase volume fill- 
ing factor of unity) is similar to what is observed for 
the dense cores of local ULIRGs dScoville et"al~l 1 1991 1: 
Downes fe Solomon! Il998t iSakamoto et al.l 119991 120081: 
Papadopoulos et ail I2010D . Furthermore, it leads to 
effective far-IR optical depths (inferred from modified 
blackbody fitting using L v cx (1 — e~ T ")B v (Td)\ see 



below for details) consistent with what is observed 
for local ULIRGs and SMGs, r = 1 at rest-frame 
A ~ 200/im for the simulations versu s at rest-frame 
A ~ 200 - 270 urn for local ULIRGs (iLisenfeld et ap 
120001 : iPapadopoulos et al.1 1201(1 iRangwala et al.1 120 lH) 
and SMGs (ILupu et all 120101 iConlev et al.l I201lh (but 
cf. iKovacs et all 12010ft . However, it is still important 
to note that the sub-resolution IS M structure is the ke y 
uncertainty in these calculations ()Younger et al.l l2009a). 
While unresolved clumpy dust can significantly affect 
the resulting SED (|Witt fc GordorJl996t lVarosi fc Dwell 
Il999ft . a more detailed treatment is beyond the scope of 
this work. The trends presented in this work should 
be qualitatively insensitive to the sub-resolution ISM 
assumption because, as explained below, the dominant 
drivers of the differences between the quiescent star for- 
mation and starburst cases are the contribution from 
stars formed pre-burst to the luminosity at the time 
of the burst and the rapid gas consumption during the 
burst, both of which do not depend on the treatment of 
sub-resolution clumpy dust. 

3. RESULTS 

Figure Q] shows the time evolution of the observed 
SCUBA 850 urn flux density (mJy), SFR (M© yr" 1 ; cal- 
culated by dividing the mass of stars formed in the last 10 
Myr of the simulation by 10 Myr), bolometric luminosity 
Lboi (Lq), initial fraction of baryonic mass that is gas f g , 
dust mass Md (M ), and dust temperaturdZH Td (K) for 
the isolated disk galaxy (left) and merger (right), where 
all quantities except f g have been normalized by divid- 
ing by their maximum values, given in the legend. When 
calculating the observed flux density we have assumed 
the simulated galaxy is at z = 2. The disk is somewhat 
unstable initially; after the disk settles, Lboi, SFR, and 
S850 decrease steadily with time. Over the 2 Gyr of the 
simulation the gas fraction decreases from 60% to 20%. 
As the gas is consumed, the SFR, and thus I/bob both 
decrease, by factors of ~ 10 and ~ 5, respectively. Ss5a 
decreases by ~ 2.5x. Md also decreases as a result of the 
decrease in gas mass, but only by ~ 40% because the de- 
crease in gas mass is partially mitigated by metal enrich- 
ment of the gas from star formation, as the metallicity 
doubles over the course of the simulation. While it may 
seem counter-intuitive that Md decreases with time, for 
a simple closed-box model assuming dust traces metals 
and c onstant yield it can be shown (jEdmunds fc Ealei 
1998) that for f g < 0.6 the dust mass increases at most 
by ~ 0.1 dex, and it decreases monotonically with f g 
for f g < 0.4. Furthermore, the preferential consumption 

11 We have calculated by fitting the modified blackbody L„ <x 
(1 - exp{-(u/u )f 3 ])B„(T d ) to the rest-frame 20 - 1000 /im SED, 
allowing all parameters to vary. Here vq is the frequency at which 
the effective optical depth t„ = 1. We have assumed the opacity 
has a power-law dependence on v in the IR, k u oc u@ . Note we 
have not used the optically thin form L v oc v& B v (T^ ) , as is almost 
always done, because we find that the first form, which does not 
assume optical thinness, provides a significantly better fit to our 
SEDs. This is because our simulated SMGs can be optically thick 
out to rest-frame > 100 jim, which i s sup ported by recent He rschel 
observations from Lupu et al. (2010) and Conlcv et al. (2011), who 
found t ~ 1 at rest-frame A ~ 200 fira. Our fitting procedure 
gives systematically higher (by as much as 20 K) than when the 
optically thin form is used, so comparisons of our dust temperatures 
to other results should take this into account. 
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Figure 1. Left: The isolated disk simulation's observed-frame integrated SCUBA 850 fim flux density (mjy, assuming z = 2; solid black) 
measured from one of the seven viewing angles, bolometric luminosity (Lq; dotted green), SFR (Mq yr ; dashed blue), dust mass Md 
(Mq; dash dot red), gas fraction f g (dash dot dot dot orange), and dust temperature (derived from SED fitting) T d (K; long dashed 
magenta) versus time (Gyr). Except for / 9 , the quantities have been normalized by dividing by their maximum values, given in the legend. 
Once the disk reaches equilibrium, Sgso, -£<bol, SFR, and f g concomitantly decrease exponentially with time as the gas is converted into 
stars. Md also decreases, but by less than a factor of 2, because the decreasing f g is offset by the increasing metallicity of the gas. Td 
decreases from ~ 60 K to ~ 50 K. Right: Same, but for the major merger simulation. All quantities are totals for the two-disk system. 
Compared to the isolated disk the time evolution is more complex: Initially 5850 1 SFR, Lboh an d Md are roughly just the sums of the 
values for the isolated progenitor disks. At first passage (t ~ 0.1 Gyr), the SFR is not elevated much beyond the baseline r ate because the 
disks a re very gas-rich and thus lack the massive stellar bar needed to efficiently remove angular momentum from the gas (IHopkins et alj 
2009b). As the two disks coalesce (t ~ 0.7 Gyr), tidal torques cause a burst of star formation, resulting in the sharp increase in the SFR 
(> 16x), Lbol (~ 7x), and S850 (< 2x) at ~ 0.7 Gyr. The increase in Lbol is much less than that boost in SFR because the luminosity 
of the stars formed during the burst is only ~ 6x the luminosity from stars already formed pre-burst. Td increases sharply from ~ 50 K 
to ~ 65 K because of the strong increase in Lb i, concurrent decrease in Md, and more compact geometry. This mitigates the increase in 
5850 that occurs from increased Lbol- The second, minor peak in Lbol, which occurs ~ 40 Myr after the peak SFR, corresponds to the 
peak AGN luminosity. Pre-coalescence, f a decreases at a rate similar to the isolated disk case. At coalescence the gas is rapidly consumed 
in the central starburst. Md decreases by a factor of 4.5, with the bulk of the decrease occurring at coalescence. 



of metal-enriched gas that occurs in our models should 
result in a lower dust mass tha n the simple closed-box 
case of lEdmunds fc Eal es (1998), which assumes perfect 
mixing. 

The behavior of the merger (Figure [T] right) is quali- 
tatively different from that of the isolated disk. Initially, 
SFR, Lbol, and Ssso are roughly equal to the sum of the 
isolated values for the two progenitor disks, because the 
disks are too gas- rich at first passage (t ~ 0.1 Gyr) for 
tidal torques to cause a strong starburst, as a significant 
stellar bar is re quired for the gas to effi ciently loose angu- 
lar momentum (|Hopkins et a l. 2009b). However, at final 
coalescence of the two disks (~ 0.7 Gyr) tidal torques in- 
duce a strong starburst, causing the SFR to increase by a 
factor of > 16. The peak of the burst is very narrow and 
significant luminosity from previously formed stars re- 
mains, so the bolometric luminosity increases by a much 
smaller amount (~ 7x) than the SFR. Moreover, as the 
gas is rapidly consumed in the starburst, Md plummets 
by a factor of 3. Along with the more compact geome- 
try, the decreased dust mass causes the SED to become 
hotter, with Td increasing from ~ 50 K to ~ 65 K. The 
increase in dust temperature during the starburst is qual- 
itatively consistent with observations, as local ULIRGs 
(i.e., merger-induced starbursts) tend to have hotter dust 
temperatures (~ 42 K) than less lumin ous (quiescent) 
galaxies (~ 35 K) (jClements et al.ll2010| ). The increased 
Td partially offsets the increase in Sgw caused by the in- 
creased luminosity. The combination of the significant 
pre-burst contribution to Lbol, the small mass of stars 
formed in the burst, and the increased Ti cause Ssso to 
increase by < 2x even though the SFR increases by > 16x 



in the burst. 

3.1. The relationship between submm flux density and 

SFR 

Figure [2] shows the observed SCUBA 850 /Ltm flux 
density in mjy versus star formation rate in units of 
Mq yr^ 1 for the isolated disk (left) and major merger 
(right) viewed from all of the 7 cameras. The submm flux 
density of the isolated disk is tightly correlated with SFR, 
increasing monotonically as SFR 4 (see best-fit curve). 
This correlation occurs because once the disk settles, Lbol 
and SFR both decrease exponentially with time. The 
dust mass also decreases, but by less than a factor of 2 
over the 2 Gyr of the simulation (see Figure [TJ. Both the 
decreased luminosity and the decreased dust mass cause 
the submm flux density to decrease. 

The case for the major merger is again qualitatively 
different. Pre-coalescence, the relationships are essen- 
tially the same as for the isolated disks. This is because 
Ss5o, SFR, Lbol, and Md at this stage are essentially just 
the sum of the two disks' isolated values, and multiplying 
all quantities by the same factor (2 for the major merger 
here) does not change the power-law index. The normal- 
ization of the relation is ~ 1.5 greater for the merging 
disks than for the isolated disk. The reason is as follows: 
An isolated disk of SFR s has Ssso = As aA , where A is 
the normalization of the SFR-S'gso relation for the iso- 
lated disk. For a non-interacting system of two identical 
disks to have total SFR equal to that of the single iso- 
lated disk, the two disks must each have SFR = 0.5s. 
Thus the total submm flux density of the system is the 
submm flux density of a single disk of SFR 0.5s, which 
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Figure 2. Integrated SCUBA 850 fim flux density (mjy) versus SFR (Mq yr x ) for the isolated disk (left) and major merger (right) 
viewed from all of the 7 different camer as. The best-fit power laws (s olid lines), linear relation ffsso = 0-5 mjy (SFR /100 Mq yr — 1 ) 
(green dashed lines; normalization f rom [Ncri et al. 2003 and Equation lAHJl , the relation for the Charv & Elbaz (2001) templates (red 



dash dot), the . Magnelli et al. (2010) relati o ns fo r their entire sample (orange dash dot dot dot) and excluding the lensed SMGs (navy 
long-dashed), the value for the IPope et al.l 1120081) composite SEP (mag enta star), and the median (cyan upward-pointing triangle) and 



mean (cyan downward-pointing triangle) from |Michalowski et al. (2010a) are also shown. For the isolated disk, the submm flux density is 
tightly correlated with both SFR and Lboli increasing essentially monotonically as SFR 4 . For the major merger, pre-coalescence (black 
open circles) the power-law index is the same as for the isolated disks because the SFR, Lbol, dust mass, and submm flux density are 
essentially two times the isolated disk values (see Figure [TJ, so only the normalization changes. During the coalescence-induced starburst 
(blue asterisks), the relationship is significantly shallower, with submm flux density scaling as SFR 01 . This is due to two effects: 1. The 
stars formed before the peak of the starburst contribute significantly to Lbol a * the starburst peak, so Lbol 9^ SFR. 2. The rapid gas 
consumption during the burst causes ALj to plummet, and the decrease in dust mass and more compact geometry cause to increase 
sharply, mitigating the increase of Ssso caused by the increased Lbol- 



we calculate using the isolated relation, multiplied by 2. 
This is 5 850 = 2A(0.5s) - 4 = l.5As 0A . Therefore the nor- 
malization of the SFR-S850 relation for the sum of two 
identical disks is 1.5 times that of the individual disk re- 
lation. This fact has important implications for the SMG 
population, which we discuss in £14.21 

On the other hand, the merger-induced burst is signif- 
icantly less effective at boosting the submm flux density. 
For a given SFR, the submm flux density is significantly 
less than for the isolated and pre-coalescence (quiescent 
star formation) cases. This is because of two reasons: 
1. The sharp decrease in dust mass and more compact 
geometry cause an increase in dust temperature, miti- 
gating the increase in 5*850 caused by increased Lbol- 2. 
The significant luminosity contributed by stars formed 
before the starburst causes Lbol to increase sub-linearly 
with SFRFl During the burst, L bo i ~ L P rc- P cak+a SFR, 
where a is the luminosity per unit star formation rate for 
an instantaneous burst. Thus Lbol is not proportional to 
SFR when L prc „ pca k is non-negligible compared to the 
luminosity of newly-formed stars, which is the case here 
because a relatively small fraction of the stellar mass is 
formed in the sharp, short-duration burst. For the burst, 
the submm flux density scales as SFR 01 (see best-fit 
curve), significantly more weakly than for quiescent star 
formation, and the ratio of submm flux density to SFR 
is significantly lower. Hence, bursts of star formation are 
significantly less effective at boosting submm flux density 
than one might naively expect. 

It is interesting to note that, during the starburst, the 
observed submm flux density can vary significantly with 

12 In principle the AGN can also cause such an effect, but for 
snapshots classified as SMGs the typical AGN contribution to the 
IR luminosity is < 10%, so the AGN is sub-dominant. 



viewing angle (e.g., for the snapshot with peak SFR, Ssso 
varies in the range ^6 — 9 mjy depending on the cam- 
era). We have confirmed that this variation is due to 
dust self-absorption: the central regions of the starburst 
can be so obscured that even the IR emission is signif- 
icantly anisotropic. As a result, the dust temperature, 
and thus submm flux density, depends on the line-of- 
sight. Though we will not explore this possibility further 
in this work, we note that differences in viewing angle 
may be enough to account for the spread of dust tem- 
peratures observed for high-z ULIRGs. In other words, 
from one viewing angle a simulated galaxy may be iden- 
tified as an SMG whereas from another viewing angle the 
same galaxy could be identified as a hot-dust-dominated 
ULIRG undetected in the submm. 

Figure [2] also shows the linear relation 5s5o = 
0.5 mjy (SFR/100 M yr" 1 ) (green d ashed lines: ob- 
tained using t he 5*8 50 — Lir relation from Neri et alj|2003 
and E quation [ATT]) , the relation for the Chary & Elbaz 
( 20011) templates (red dash dot), the Ma gnelli et al. 
(2010) relations for their entire sample (orange dash 
dot dot dot) and excluding t he lensed SMGs (navy 
long-dashed), the value for the IPope et al.l ([2008) com- 
posite SED (magenta star), and the median (cyan 
upward-pointing triang le) and mean (cyan down ward- 
pointing triangle) fro m [Michalowski et al ] (l2010af) . The 
iCharv fc Elba^ (|200l and lPope et al l ([20081 ) values were 
obtained by redshifting the templates to z — 2 and con- 
verting Lir of each t empl ate to SFR using EquationlATTl 
We u sed Equation IA11I to convert the Magnelli et al. 
(|2010D relations from L m to SFR. 

The typical va l ues fro m IPope et all (j2008T ) and 
I Michalowski et all (|2010af ) are consistent with the data 
from our major merger simulation. As explained above, 
the relations we find are much shallower than linear, 
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so the iNeri et all (|2003f) relation differs significantly 
from our relations f or both the qu i escent and star- 
burst modes. The iCharv fc Elbazl (|2001ft templates 
are also very discrepant because high-redshift ULIRG 
SEDs are often better fit by local template s appropri- 
ate f o r less luminous, co l der galaxies (e.g . . iPope et al.1 
2001 iDannerbauer et all [Mot iRex et all I2010H . The 
Mag nelli et al.l (|2010f) relations agree better with our sim- 
ulations: for the full sample, S^o oc SFR 48 . This is a 
slightly steeper relation than what we find for quiescent 
disks and significantly steeper than that for starbursts. 
Wh en the six lensed SMGs are removed from their sam- 
ple, IMagnelli et all find S 850 oc SFR 29 . The power- 
law index of this relation is less than that for our quies- 
cent disks but greater than that for our starbursts. The 
lensed SMGs tend to be intrinsically fainter and thus less 
likely to be strong starbursts than the non-lensed pop- 
ulation, so it is reasonable that inclusion of the lensed 
SMGs le ads to a steeper relation. While it is interesting 
that the IMagnelli et al.l relation for the unlensed SMGS 
is crudely consistent with what we expect for a mixed 
population, one should not over-interpret this compari- 
son. As we will discuss below, Sgso cannot be determined 
solely from SFR because the dust mass plays a signifi- 
cant role also. However, a robust conclusion that should 
be drawn from Figure [5] is that, in the simulations, the 
starburst mode is less efficient at boosting submm flux 
than the quiescent mode. 

3.2. Dependence of submm flux density on SFR, Lbob 

and Md 

For a given SFR, galaxies of different masses tend to 
also have different dust masses; thus the normalization 
of the S850-SFR relation varies for different mass sim- 
ulations though the scalings are similar. As a result, 
one cannot calculate the submm flux density given only 
the SFR, but it is possible to parameterize the submm 
flux density as a function of SFR and dust mass. Since 
much of the discrepancy in the .Ssso-SFR relations for 
quiescent star formation and starbursts is caused by the 
rapid decrease in dust mass during the starburst, we ex- 
pect that including dust mass in our parameterization 
will eliminate much of the difference between quiescent 
and starburst star formation modes. We have analyzed 
the fu ll set of simulations fr om our SMG number counts 
work (jHavward et al.l I2011L C. Hayward et al. 201 f, in 
preparation), which includes a range of progenitor disk 
baryonic masses (~ 3.5 x 10 10 — 4x 10 11 Mq), mass ratios 
(~ 0.1 — 1), and initial gas fractions (0.6 — 0.8), fitting 
the submm flux density as a power law in both SFR and 
dust mass. Both the quiescent and starburst phases are 
included. Perhaps surprisingly, the following relations 
(for simulated galaxies placed at redshift z = 2) hold to 
within ~ 0.1 dex for all but a few outliers over the range 
0.5 mJy < Ssso < 15 mJy: 



0.58 




S 85 o =0.65 mJy 
Si.i=0.30 mJy 



where Ss5o and S%,i are the fluxes in the SCUBA 850 fim 
and AzTEC 1.1 mm bands, respectively. 



Note that the SFR exponent in these relations is sim- 
ilar to that for the ffssp -SFR relation for quiescent SF 
(~ 0.4). As explained in £13 . H the sharp decrease in dust 
mass during the starburst is one of the main reasons the 
"S850-SFR relation is much shallow for starbursts than 
for quiescent SF. Adding Md as a parameter effectively 
decouples this effect; in other words, for fixed dust mass 
the S850-SFR relations for the two SF modes are much 
more similar than when the evolution of the dust mass 
is taken into account. Remaining differences caused by 
the contribution from stars formed pre-burst to the burst 
luminosity, AGN contribution, geometry, and other fac- 
tors prevent Eq. {1} from recovering Sgso exactly, but 
the small scatter suggests that these factors are subdom- 
inant. 

We can also fit the submm flux density as a function 
of L b oi and M d : 



5 850 = 0.40 mJy 



St. i=0.20 mJy 



T \ 0.52 / „ „ \ 0.60 

^bol \ / Md 



10 12 L C . 

Lho\ 
10 12 L, 



V 10 8 M Q 

( M d 
V10 8 M Q 



(2) 



0.63 



These fitting functions are accurate to within ~ 0.15 dex. 
Replacing Lboi with the bolometric IR luminosity .Lir 
yields a similar result. Figure |3] shows how well these 
fitting functions reproduce the submm flux density of 
our simulated galaxies. It is discussed in more detail in 

EH 

3.2.1. Relations for an optically thin modified blackbody 

It is instructive to compare the above relations to those 
for a single-temperature mass of dust transparent to its 
own emission, the model which is implicit in the stan- 
dard method of fitting a modified blackbody to the IR 
SED. For a mass of dust in thermal equilibrium with 
temperature Td we can express the submm flux density 
as a function of dust bolometric luminosity Ld and dust 
mass Md or SFR and Md- Assuming z — 2 and far-IR 
spectral index /3 = 2, the relations are (see the Appendix 
for a derivation) 



S 85 o = 1.4 mJy 



SFR 



100 M Q yr" 1 



1/6 



/ M d 
V 10 8 M Q 



5/6 



and 



£ 8 5o = 1-4 mJy 



10 12 L 



1/6 



( M d 



V1O 8 M 



5/6 



(3) 



(4) 



These should be compared to Equations JT]) and @, re- 
spectively, by assuming L^ \ ~ Ld (i.e., the luminosity 
emitted by stars and AGN is completely absorbed by 
dust), which is a reasonable approximation for snapshots 
classified as SMGs (5*850 > 3 mJy). The submm flux de- 
pends only weakly on redshift for the redshift range of 
interest (1 < z < 5), scaling as (1 + zf^D^ 2 , where 
Da is the angular diameter distance at redshift z (see 
Equation lA"9j) . 



3.2.2. Comparison of the relations to the full radiative 

transfer 
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Figure 3. Left: Logarithm of the ratio of the submm flux density calculated using one of the analytic forms (/Sgso.analytic) to the submm 
flux density calculated through the full RT (Sgso.Rr) vs - S'sso.RT for all time snapshots of our simulated galaxies. Black circles show the 
ratio when Ssso. analytic i s calculated from the SFR and of our simulated galaxies using Equation {T}. Blue triangles show the ratio 
when the optically thin modified blackbody model (Equation [3j is used. The line Ssso, analytic = -S'sso.RT is shown to guide the eye. Flight: 
Same, except now the black circles show the values when Equation 10 is used to calculate the submm flux density from Lbol an d ^d °f 
the simulation snapshots, and the blue triangles show the values when Equation Q is used, assuming Lbol ~ ^d- In both cases the simple 
optically thin modified blackbody model overpredicts the submm flux density by > 0.3 dex for the simulated SMGs, and the overprediction 
is worse for lower Ssso- 



Figure [3] shows the logarithm of the ratio of the 
submm flux density calculated using the above equa- 
tions (Ss5a, analytic) to the submm flux density calculated 
through the RT (S'sso.RT) versus Ssso, rt- The left panel 
shows the results when Equations (HJ black circles) and 
([3J blue triangles) are used to calculate Sgso from SFR 
and Md- The right panel shows the results when Lbol 
and Md are used instead, with the black circles corre- 
sponding to Equation ([2|) and the blue triangles Equation 
|4]). The fitting functions derived from the simulations 
are able to reproduce S^so from SFR (Lbol) and Md to 
within ~ 0.1 (0.15) dex. The simple modified blackbody 
model tends to overpredict the submm flux density by 
> 0.3 dex; the typical over-prediction is ~ 0.5 dex (a 
factor of 3). (Note that the corresponding uncertainties 
in SFR, Lbol, and Md for fixed observed submm flux den- 
sity would be less because the submm flux density scales 
with these quantities sublinearly.) Furthermore, the er- 
ror in the prediction correlates with SFR (Lbol) and dust 
mass because of the differences in the power-law indices 
for the fitting functions and the modified blackbody re- 
lations. 

The optically thin modified blackbody model fails for 
multiple reasons. From SED fitting we find that the sim- 
ulated galaxies can have effective optical depth r > 1 out 
to rest-frame ~ 200 /Ltm. For fixed dust temperature, the 
optically thin assumption will result in an overestimate of 
the luminosity density at wavelengths for which r > 1 be- 
cause (1 - exp[-(t//t/ ) /3 ]) < {v/vqY for all v > 0. Thus 
Equation (|A3|) will overestimate Ld for fixed Td and Ma- 
li Ld and Md are fixed, the dust temperature will be 
underestimated when optical thinness is assumed, and, 
therefore, the submm flux density will be overestimated. 
Also, Ld ~ Lbol is less accurate an approximation for the 
faintest sources tha n for the brightest b ecause L^/Lboi 
increases with Lbol (|Jonsson et al J 12006). If Ld < Lbol, 
the assumption that Ld ~ Lbol will overestimate Ld and 
thus overestimate the submm flux density; this may ex- 
plain why the overprediction is worse for lower S^so- Fi- 



nally, we have seen above that the assumption Lbol oc 
SFR is invalid during the burst. 

4. DISCUSSION 

We have demonstrated that the submm flux density of 
a galaxy scales differently with SFR for quiescent star 
formation and starbursts. The ratio of sub-mm flux 
density to SFR is significantly less for merger-induced 
bursts than for quiescent star formation. This is because 
of the rapid decrease in dust mass and more compact 
geometry during the starburst, which causes the SED 
to become hotter, and the significant contribution from 
stars formed pre-burst to the luminosity during the burst, 
which makes the luminosity increase by a much smaller 
factor than the SFR. As a result, merger- induced star- 
bursts are less efficient at boosting submm flux density 
than one might naively expect. 

Our results have a number of important implications; 
we discuss these now. 

4.1. Predicting submm flux densities from models 

One implication of this study is that, at a fixed red- 
shift, the galaxies with highest submm flux density are 
not necessarily those with the highest bolometric lumi- 
nosities or SFRs. Thus theoretical models, be they sim- 
ulations or semi-analytic, must explicitly calculate the 
submm flux density of their simulated galaxies in order 
to select which are SMGs as opposed to simply select- 
ing the most rapidly star-forming or most luminous ob- 
jects. However, the computational expense required to 
self-consistently calculate the submm flux density limits 
this approach to idealized, non-cosmological simulations 
(as done here), individual galaxies excised from cosmo- 
logical zoom-in simulations, or semi-analytic models in 
which various simplifying assumptions must be made. 
As an alternative to performing RT, cosmological sim- 
ulations and semi-analytic models can use the relations 
among submm flux density, SFR or bolometric luminos- 
ity, and dust mass presented herein (Equations [T] and [2]) 
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to assign submm flux density to simulated galaxies. Ad- 
ditionally, observers can use the relations to estimate the 
instantaneous SFR given submm flux density and dust 
mass (obtained from fitting the IR SED using the full 
modified blackbody as we have done or by measuring 
the gas mass and assuming a dust-to-gas ratio; ignoring 
uncertainties on submm flux density, a dust mass accu- 
rate to a factor of 2 gives SFR accurate to a factor of 
3). 

4.2. Heterogeneity of the SMG population 

These results also imply that the SMG population is 
heterogeneous. We have seen that it is possible for a 
quiescently star-forming disk to have submm flux density 
equal to that of a merger with much higher SFR (Figure 
[2]). Furthermore, since the scaling of submm flux density 
with SFR is sublinear, adding two equal disks (and thus 
doubling the dust mass and SFR of the system) increases 
the submm flux density more than simply boosting the 
SFR by a factor of 2; Equation (TT]) shows that SFR would 
have to be boosted by 5x to achieve a 2x boost in submm 
flux density if dust mass is kept constant. When the 
sharp increase in dust temperature during the starburst 
and the narrowness of the burst are accounted for the 
effect becomes even stronger: in Figure Q] we see that a 
> 16x increase in SFR gives a < 2x increase in submm 
flux density. 

The single-dish submm telescopes used for wide-field 
surveys of SMGs have beam sizes > 15"(> 130 kpc at 
z = 2). Thus, during a merger the two progenitor galax- 
ies will spend a considerable amount of time within the 
area of the beam. From the above arguments, we see 
that this is a very efficient way to create an SMG, but 
this contribution has been relatively unappreciated. We 
argue that the SMG population attributable to merg- 
ers is bimodal: some are merger-induced starbursts and 
some are two (or more) infalling disks (normal galaxies 
that are not yet interacting strongly) blended into one 
submm source ("galaxy-pair SMGs")E3 Note also that 
not only major me rgers but also fayorably oriented minor 
mergers (see, e.g., ICox et al.l[2008tlHopkins et al.H2009bl) 
can contribute to the SMG population. 

Furthermore, the most massive, rarest 'isolated disks', 
and even small groups, may also contribute; we expect 
this contribution to be subdominant because SMGs are 
on the exponential tail of the mass function, but we defer 
a precise determination of their contribution to future 
work. Though both the merger-induced starburst and 
galaxy-pair populations are mergers, only in the former 
is the star formation merger-driven (and only partially, 
as the baseline star formation that would occur in the 
disks even if they were not interacting is significant). 

13 1 Wang et alj 120111 ) recently presented high-resolution submm 
continuum images of two SMGs which were previously identified 
as single sources but are resolved as 2 or 3 distinct submm sources 
in their images. The sources are at significantly different redshifts 
and thus physically unrelated. Our galaxy-pair SMGs are also two 
sources blended into one, but they are distinct from the type of 
SMGs Wang et al. observed because they are merging and thus 
physically connected. Both types of blended sources are poten- 
tially important SMG subpopulations that complicate our under- 
standing of SMGs, and it is crucial to understand the relative con- 
tributions of merger-induced starbursts, galaxy pairs/infall-stage 
mergers, quiescently star-forming disks, and physically unrelated, 
blended sources. 



Given that the physically meaningful property of local 
ULIRGs is not that their IR luminosities are > 10 12 L Q 
but that they are powered by merger-driven star forma- 
tion and AGN, only the merger-driven starburst category 
of SMGs should be considered physically analogous to lo- 
cal ULIRGs. 

The observational signatures and physical implications 
of this bimodality will be discussed in future work. 
Here we simply note that the galaxy-pair contribution 
is suppor ted observationally by the frequency of multi- 
ple radio (llvison et alJl2002L [2007t iChapman et al.ll2005l : 
lYouneer et al.l l2009d) . 24 /urn (|Pope et all l200ffl . and 
350 /Ltm (|Kovacs et al.1 120101 ) counterparts to SMGs; 
by CO interferometry showi ng that a large frac t ion of 
SMGs are resolved b i naries dTacconi et all 120061 120081: 
iBothwell et all [20101 lEngel et all I2O1O0 : and by the 
SMGs that have mo rphologies that do not resemble 
merger remnants (e.g.,jB othwcll et al. 2010; Carilli ~et al.1 
EMU IRicciardelli eTafll2lTi7il : iTargett et alllMTll) 

4.3. SMG masses 

The masses of SMGs are hotly debated, with differ- 
ent authors finding masses discrepant b y ~ 6x for the 
same SMGs (jMichalowski et al.l l2010al lbl: Mainline et al.1 
I2010T) . Accurate masses are important in order 
to test potential evolutionary relatio nships among 
SMGs and other galaxy cl as ses (e.g.. iBrodwin et al.l 
2008t iBussmann et al.ll2009aHb1 : iNaravanan et al.l l2010b: 
Rothberg fc Fischerll2010[) " and to check that number den- 
sities of SMGs are consistent with observed stellar mass 
functions. Stellar mass determinations from SED fitting 
are limited by uncertainties in stellar evolution tracks, 
the initial mass function, star formation histories, dust 
attenuation, and AGN contamination. Since our models 
use star formation histories, attenuation (from the ge- 
ometry of stars, AGN, and dust), and AGN components 
that originate directly from the hydrodynamic simula- 
tions instead of the standard assumptions — e.g., instan- 
taneous burst or exponential star formation histories , 
the Calzetti a ttenuation law (|Calzetti et al.l 11994 120001: 
lCalzettilll997|) — we can provide constraints on SMG stel- 
lar masses that are complementary to those derived from 
SED fitting. 

Given the inefficiency of bursts at boosting submm flux 
density that we have demonstrated above, SMGs must be 
very massive, because smaller galaxies undergoing even 
very strong bursts cannot make SMGs. Our models re- 
quire M* > 6 x 10 10 M o to reach S S50 > 3 mJy, and 
typical masses are higher. The area of the Ss5o — M+ 
plot spanned by our models agrees well with the obser - 
vationally derive d values of | Mic halowski "et~aLl ((2010a). 
About half of the Mainline et al.l (|2010f ) values lie in the 
area spanned by our models, whereas the other half have 
lower masses. However, the single-component star for- 
mation histories assumed by Hainline et al. may cause 
the stellar masses to be underestimated by ~ 2x, which 
would resolve much of the discrepancy. A detailed com- 
parison of the mass estimates will be presented in M. 
Michalowski et al. (2011, submitted). 



14 Note howe ver, that the IBothwell et"aU ^20101) and 
ICarilli et al. ( 2010'j) objects that resemble disk galaxies may in fact 
be the molecular disks that re-form rapidly after a ga s-rich merger 
IjNaravanan et al.|[2008bl : [Rob ertson & Bullock 2008). 
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4.4. SMG duty cycles 

Understanding the duty cycle of SMGs is important 
for predicting submm counts from models, quantifying 
the contribution of SMGs to stellar mass buildup, and 
interpreting star formation efficiencies of SMGs. Since 
the submm flux density depends on luminosity heating 
the dust, dust mass, and geometry, the submm duty cy- 
cle depends on the same factors. As we have seen, the 
starburst induced at merger coalescence causes a sharp 
peak in SFR, Lboi, and submm flux density. However, 
the duty cycle is limited because of the sharp cutoff in 
SFR, and thus drop in Lboi, after the burst and the sig- 
nificant drop in dust mass that occurs as highly enriched 
gas is consumed in the burst. 

Since the submm flux density depends more strongly 
on dust mass than either SFR or Lboi (see Equations 
Q]and[2]), it i s more efficient if one can keep more dust 
around at the expense of lower star formation rate. The 
quiescent star formation mode does exactly this. As a 
result, the galaxy-pair phase (discussed in fl4.2j) adds sig- 
nificantly to the submm duty cycle. Figure Q] shows that 
the galaxy-pair phase has a longer submm duty cycle 
than the burst, though the lack of smooth accretion in 
our simulations — and thus need for starting with very 
high initial gas fractions — complicates a precise determi- 
nation of the relative flux densities and duty cycles of the 
galaxy-pair and starburst phases. Regardless, it is clear 
that the galaxy-pair phase increases the SMG duty cy- 
cle significantly, alleviating some of the tension between 
the submm counts estimated from high-redshift major 
me rger rates and shor t (~ 100 Myr) starburst duty cycles 
by iDave et"al~l (|2010l ) and the observed submm counts. 
Inclusion of both SMG populations is crucial to match 
the observed SMG number coun ts without resorting t o a 
top-heavy initial mass function (Havward ct al. 20 111 C. 
Hayward et al. 2011, in preparation). 

4.5. Implications for IR SED fitting 

One reason the single-temperature, optically thin mod- 
ified blackbody fails is that the effective optical depth 
of our simulated SMGs can be greater than 1 out to 
rest-frame ~ 200 /im. This is consistent with th e ef- 
fective optica l depth s derived by iLupu et al.l (|2010l ) and 
iConlev et al.1 (|2011l ) when they fit a general modified 
blackbody (i.e., not assuming optical thinness) to the IR 
SEDs of their SMGs. By assuming optical thinness and 
only fitting longward of the FIR peak, one overestimates 
the lumino sity density at wavelengths for which r > 1 
(see §3.2.2[) P°1 Consequently, the assumption of optical 
thinness yields a colder dust temperature (by as much as 
~ 20 K) than if optical thinness is assumed. 

In the pve-Herschel era, often flux densities in only a 
few FIR bands and at 1.4 GHz were available for large 
samples of SMGs. As a result of the limited number of 
data points, models more complex than the optically thin 
modified blackbody (e.g., the full modified blackbody 
used here or models assuming a distribution of temper- 
atures) cou ld not provide a be tter description of the IR 
SEDs (e.g. JKovacs et al.ll2006l ). Recently, some authors 



15 This may explain why IPope et al.l J2006) found that the 
submm flux density tends to overpredict LrR, and the overpredic- 
tion is worse for SMGs at low redshift where the SED is sampled 
further from the IR peak. 



have used mode ls incorporating a d i stribution of tem- 
peratu res (e.g., IKovacs et al.l 120101 : iMic halowsk i et all 
l2010al) . fin ding that such mod els better described the 
IR SE Ds (IKovacs et al.1 12010ft Now that Hersche l 
PACS (iDannerbauer et al.l 120101: [Magnelli et al.l l2010f > 
and SPIRE |Chapman et al.1 12010 ) have provided data 
over the entire FIR SED for large samples of SMGs it is 
possible to perform more sophisticated fitting for many 
SMGs. Given the physical inferences that are drawn from 
effective dust temperatures obtained via FIR SED fitting, 
it is important to have as robust a method as possible 
and to take full advantage of the available data. We will 
present such a method in future work. 

4.6. Limitations of our model 

At this point we find it instructive to define the lim- 
itations of this work so that the results can be placed 
in an appropriate context and future experiments can be 
designed for maximal impact. One of the primary limi- 
tations, both for the radiative transfer and the hydrody- 
namics, is the treatment of the sub-resolution interstellar 
medium, especially — because we focus upon the submm 
flux — the structure, distribution, and composition of the 
dust which dominates emission at submm wavelengths. 
In fact, the differences between the model employed in 
this work and that we used in N10 were motivated by a 
desire for a simpler treatment of the sub-resolution ISM 
(see §2.2.1[) . In N10, stars with age < 10 Myr dominated 
the submm flux, so the submm flux was closely tied to 
the SFR (see Fig. 1 of N10). Our simplified assumptions 
(no sub-resolution PDR model, uniform ISM density on 
scales below the SPH smoothing length) result in submm 
flux that is tied more directly to the bolometric luminos- 
ity than the instantaneous star formation rate because 
stars with age > 10 Myr contribute significantly to the 
bolometric luminosity and, because they are still deeply 
embedded in dust, the submm flux. 

While our model has the advantages of simplicity and 
strict physical consistency (because the obscuration orig- 
inates purely from the hydrodynamic simulations rather 
than from sub-resolution PDRs), we still must make 
an assumption about the sub-resolution structure of 
the ISM. By assuming uniform density on scales be- 
low the smoothing length we only include dumpiness 
that arises from the hydrodynamic simulations, so this 
assumption may be considered conservative. The as- 
sumption is also simplistic, of course, because the real 
ISM has significant structure on scales < 100 pc. How- 
ever, proper treatment of radiative transfer through a 
clumpy medi um is a significant area of research in and 
of itself (e.g.. iHobson fc Padmanl [l993t iWitt fc Gordon! 
119961 : IVarosi fc Dweklll999D and thus beyond the scope 
of this work. A study of the effects of sub-resolution 
dust dumpiness on galaxy SEDs and efforts to devise 
a better treatment of sub-resolution dust dumpiness in 
Sunrise are underway. 

We also caution that modifications to the simple star 
formation prescription and ISM treatment in the hy- 
drodynamical simulations themselves could change the 
amplitude and dura tion of starbursts (|Cox et al.l l2006b: 
iSpringel et al.l [20051) . However, while changing the SF 
prescription or ISM treatment could change the rela- 
tive contribution of quiescent and starburst star for- 
mation modes to the star formation history of a given 
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merger, this alone should not change the differences be- 
tween quiescent star formation and starbursts which lead 
to the significantly different relationships between SFR 
and submm flux. Changing the SF or ISM prescription 
could substantially alter the spatial distribution of dust 
and stars and thereby modify the detailed relations be- 
tween SFR/Lboi, Mdi and Ssso- However, because ge- 
ometry is relatively unimportant in setting the relations, 
this uncertainty should have a relatively limited influ- 
ence on our results. Furthermore, changing the feedback 
implementation can alter the evolution of the IR SEDs 
(jChakrabarti et all 120071) . On going and future studies, 
with much higher resolution and more advanced track- 
ing of the clumpy ISM (e.g., iHopkins et al|[20TTI ). will 
improve the predictive power of our models. 

Furthermore, we stress that the simulations presented 
here are not cosmological. We adopt this approach be- 
cause it enables us to achieve the high resolution neces- 
sary to perform radiative transfer in order to accurately 
calculate the submm flux density; to survey the param- 
eter space of progenitor masses, mass ratios, and orbits; 
and to avoid uncertainties in modeling realistic galaxy 
populations in a cosmological environment. The primary 
drawback of this approach for our present purposes is the 
lack of gas accretion, which cosmological hydrodynamic 
simulations show is a significant driver of star formation 
for the high -redshift, massive galaxies with which we ar e 
concerned (jKeres et all 120051 120091 : iDekel et all 120091 ) . 
Gas accretion can continually replenish the gas in the 
galaxy, maintaining relatively high gas fracti ons and rel- 
ative ly constant star formation histories (e.g. JDave et all 
120101) . 

Inclusion of cosmological gas accretion would alter the 
time evolution presented in Figures [TJ but it would not 
significantly alter the differences between the quiescent 
and merger-induced burst modes of star formation. This 
is because mergers would still superimpose a strong burst 
of star formation and sharp decrease in gas mass over the 
baseline evolution. Furthermore, unless smooth accre- 
tion significantly affects the geometry of stars and dust 
in the simulated galaxies, it will not have a significant 
effect on the relationship between submm flux, Lboi, and 
dust mass (the relationship between submm flux, SFR, 
and dust mass may be more affected because the relation 
between SFR and Lboi may be changed significantly). 
Thus inclusion of cosmological gas accretion should not 
qualitatively alter our results. 

4.7. Connections to previous work 

In previous work, we developed a model relat- 
ing the evolution of ga laxies, starbursts, and quasars 
(IHopkins et all l2006al fbl l2008al l2009dl : iSomerville et all 
2008). A principal conclusion from these analyses is 
that while starbursts driven by gas-rich mergers can ac- 
count for many instances of unusual activity in galaxies, 
they provide only a minor contri bution to the star for - 
mation history of the Univer se (Hopki ns et al.l l2006cf ). 
Indeed, IHopkins et all (|2010l ) emphasize that much of 
the star formation during galaxy interactions occurs 
in the "quiescent" mode and should not be counted 
as part of a merger-induced starburst. This is sup- 
ported by the decompo sition of the light p rofiles of 
nearby o ngoing mergers (IHopkins "elall l2008bh and lo- 
cal cusp (|Hopkins et al.ll2009aD and core ( Hopkin s et all 



l2009d ) ellipticals. All of t hese objects exhibit evi- 
dence of "excess" centr al light (jRothberg fe Josephll2004l ; 
iKormendv et al.l I2009D. indicativ e of relic starbursts 
(jMihos fc Hernquistl II99I I1996D caused by merger- 
driven inflows of gas ([Barnes fc Hernquistl 1 1 9911 [1996). 
The integrated mass in these components agrees well 
with estim ates of the cosmic history o f merger-induced 
starbursts ()Hopkins fc Hern quist 20l3). 

The results presented herein extend these conclu- 
sions to high-redshift phenomena. Critically, we find 
that "quiescent" star formation during galaxy interac- 
tions, i.e., star formation which occurs during the in- 
fall/pair stage, is a key element in understanding the 
brightest submm sources, especially their number counts 
and duty cycles, and connecting them to other high- 
redshi f t populations including quasars (Hopkin s et all 
12008a!: iNaravanan et al.ll2008bl) a nd compact spheroidal 
galaxies (jWuvts et alll2009U2010l ). Just as nearby LIRGs 
represent a heterogeneous collection of merging and iso- 
lated systems, it is natural to suggest that the population 
of high-redshift submm galaxies is heterogeneous, as we 
have argued here. 

5. CONCLUSIONS 

We have combined high-resolution 3-D hydrodynamic 
simulations of high-redshift isolated and merging disk 
galaxies and 3-D Monte Carlo dust radiative transfer 
calculations to study the sub-millimeter galaxy selection, 
focusing on the relationships among submm flux den- 
sity, star formation rate, bolometric luminosity, and dust 
mass. Our main conclusions are the following: 

1. The relationship between SFR and submm flux 
density differs significantly for quiescent and star- 
burst star formation modes. Starbursts produce 
significantly less submm flux density for a given 
SFR, and the scaling between submm flux density 
and SFR is significantly weaker for bursts than for 
quiescent star formation. Bursts are a very inef- 
ficient way to boost submm flux density (e.g., a 
starburst that increases SFR by > 16x increases 
submm flux density by < 2x). Another conse- 
quence is that the galaxies with highest submm flux 
density are not necessarily those with highest SFR 
or bolometric or infrared luminosity. 

2. The submm flux density of our simulations can be 
parameterized as a power law in SFR and dust 
mass (Lboi and dust mass) to within ~ 0.1(0.15) 
dex. The scaling derived from the commonly used 
optically thin modified blackbody model systemat- 
ically overpredicts the submm flux density by > 2x 
because numerous assumptions of the model (opti- 
cal thinness in the FIR, Lir oc SFR, Lboi ~ Lir) do 
not hold. The fitting functions we provide (Equa- 
tions [T] and [5]) should be useful for calculating the 
flux density in semi-analytical models and cosmo- 
logical simulations when full radiative transfer can- 
not be performed and for interpreting observations. 

3. Mergers create SMGs through another mecha- 
nism besides the strong starburst induced at 
coalescence — they cause the two infalling disks to 
be observed as one submm source because both 
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disks will be within the large (~ 15" , or 130 kpc at 
z = 2) beam of the single-dish submm telescopes 
used to identify SMGs during much of the infall 
stage. For major mergers, this effect boosts the 
submm flux density by 2x. To achieve the same 
boost in submm flux density one would have to 
boost the SFR of a quiescent disk by ~ 6x or in- 
duce a starburst that boosts the SFR by > 16x. 
This implies that the SMG population is hetero- 
geneous: it is composed of both late-stage major 
mergers and two (or more) infalling disks observed 
as a single submm source ("galaxy-pair SMGs"). 
The largest quiescently star-forming galaxies may 
also contribute. Thus, unlike local ULIRGs, SMGs 
are a mix of quiescent and starburst sources. 

4. SMGs must be very massive: to reach 5850 > 3 
mJy, stellar mass of at least 6 x 1O 1O M is required, 
and typical values are higher. 

5. The submm duty cycles of our simulated galaxies 
are a factor of a few longer than what one would 
expect if all SMGs were merger-driven bursts be- 
cause the relatively gentle decline in SFR, Lboi, and 
dust mass during the galaxy-pair phase results in 
a longer duty cycle for the galaxy-pair phase than 
for the starburst. The duty cycle of the latter is 
limited because the peak in luminosity is narrow 
and the dust temperature increases sharply during 
the burst. 

6. Fitting the SEDs of SMGs with an optically thin 
modified blackbody tends to yield significantly 
lower dust temperatures than when the full opac- 
ity term is used because the effective optical depths 
can be ~ 1 out to rest-frame ~ 200 /jm, both for 
our simulated SMGs and observed SMGs. There- 
fore, one should be cautious when interpreting ef- 
fective dust temperatures derived via fitting an op- 
tically thin modified blackbody to the FIR SEE), 
especially when comparing SMGs to galaxies for 
which optical thinness in the IR may be a reason- 
able approximation. 

Future work will include predictions of submm num- 
ber counts from our model, an investigation of the ob- 
servational signatures and physical implications of the 
proposed SMG bimodality, and an improved method for 
fitting IR SEDs of galaxies. 
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APPENDIX 

DERIVATION OF THE RELATIONS GIVEN IN §3.2.1 

Here we derive the relations for submm flux density as a function of dust bolometric luminosity Ld and dust mass 
M d (Equation HI) and SFR and Md (EquationBl) for an optically thin modified blackbody. One can mode l galaxy SEDs 
with more complex models fe.g.. lDale et~alll2001l : iDale fe Heloull2002l : iChakrabarti fe McKedl2005l [2008HKovacs et al.1 
2010), but for the sake of simplicity and because the optically thin modified blackbody is commonly used for SED 
fitting we only consider an optically thin modified blackbody here. Consider a mass Md of dust with temperature T d . 
Assuming the dust is optically thin at rest- frame frequency v r , the luminosity density emitted by the dust at that 
frequency is 

L„ r =4nK Vr M d B Vr (T d ), (Al) 

where K Ur is the dust opacity (m 2 kg^ 1 ) at rest-frame frequency v r and B U7 ,(Td) is the Planck function. We assume a 
power-law opacity in the IR, 

K Vr = K , (A2) 

where kq is the opacity at frequency Vq. Integrating Equation (IA1|) over v gives the total dust luminosity, 

U = F(4 + C (4 + 0) V M dK0 (*< Y , (A3) 



c 



Kvq ) 



where T and ( are Riemann functions, h is the Planck constant, c is the speed of light, and k is the Boltzmann constant. 
Solving for the effective dust temperature yields 



h 

Td =k 



2,.P 



Ld&vP 



l/(4+/3) 



(A4) 



T(4 + /3)C(4 + l3)&TTK hM d 
If we place the mass of dust at redshift z, the flux density at observed- frame frequency v is 

S Vo (T d ) = (l + z)^ (A5) 

= {l + z) ^M d B„ r {T d ) (A6) 

= {1 + z f-^f^Y B {Td)t (A7) 

where we have related angular diameter distance Da and luminosity distance Dl using Dl = (1 + z) 2 Da- In the 
Rayleigh-Jeans limit, B„(T) = 2kv 2 T/c 2 , so 

S Vo {T d ) = {l + zf- 1 ^ ( -V ' vlMdTd. (A8) 



c 2 D\ \u 

By substituting Td from Equation (|A4|) into Equation (|A8j) we find 

/ \ 9 / 2 \ 1 /( 4 +^) 

2hK ..2 ( V o \ I VqC \ ^ i /(4+/ 3) „(3+/3)/(4+/3) 



s »°- c 2 d 2 /°{vJ [r(A+m^+m^oh) {1 + z} Ld Md ■ (A9) 

For the I Weingartner fc Dra ine (2001) Ry = 3.1 Milky Way dust model, which we use in our RT calculations, /3 w 2 
and the 850 ^m opacity is k 8 so = 0.050 m 2 kg -1 , consistent with the value Uames et al.1 (|2002l ) derived from submm 
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observations of local galaxies, 0.07 ± 0.02 m 2 kg -1 , and with the results of iDunne et al.l ()2003l ). Thus the observed 
850 fim flux density is 

*- = " <> + *> (t^)" 2 (to^)'" {mk)"'- < A10 > 

The iKennicuttl (|1998bf ) SFR-Lir calibration converted to a iKroupal (|2001l ) initial mass function is 

Lm « Lbol « 9 x 1O 9 L (SFR/M yr" 1 ). (All) 

This conversion assumes all starlight is absorbed by dust and the contribution from AGN and old stars is negligible; 
as discussed above, these assumptions are all violated at some level. (If these assumption were true, the power-law 
indices in Equations ([1]) and ([5} would be identical.) However, since the above calibration is ubiquitously applied, we 
will give the relation that results when we use it; the relation is 

Assuming Q m = 0.27, = 0.73, and h = 0.7, the angular diameter distance at z = 2 is 1.77 Gpc H Wright! 12006ft . 
Thus for z = 2 we recover Equation 



.... / UO 8 M 0/ 



^o = 1.4mJy(-^V / 7-^-) . (A13) 



This should be compared to Equation (j2J. In terms of SFR, we get Equation (j^J), 

S S5Q = 1.4 mJy ( SFR \ 1/6 / _m\ 5/6 ( } 

850 y \100 M Q yr" 1 / \10 S M Q J V ' 

This should be compared to Equation ([!])■ 

Even if the underlying power- law index of the dust opacity curve is ft = 2, for a dist ribution of dust tem- 
peratures a single-temperat ure modified blackbody with = 1.5 may better fit the SED (jDunne fc Ealesl I200H : 
iChakrabarti fc Mc Kce 2008ft . Note also that the fitted dust temperature and are degenerate, and t he fitted val- 



ues ca n depend sensitively on both noise in the data and temperature variations along the line-of-sight (Shetty et al 



2009 ap. S ince = 1.5 is often assumed when fitting SEDs and determining dust masses of SMGs (e.g.. lKovacs et al 
20061 boidt ICoppin et alll2008t iChapman etaHl2010ft . so we will provide the relations for = 1.5 also. They are: 



n \ — 2 / r \ 0.18 / „ , T \ 0.82 
1->A \ I L d \ f Md 



^ = 1.9^(1 + ^1^ [j^J (A15) 



and 



,s-,„ = 1-9 mJy (1 + z f- 5 ( 2 ( SFR V' 18 ( _M±_ ) . (A16) 

Sta = 1 . 0m3y{ j±-y i, (^±-y*\ ( Ai7) 



For z 



and 



SFR \ ul ° / M, 



g850 = 1.0 mJ y^ iQO — yr _ x j ^-j . (A18) 

The equations can be rescaled to different values of kq using S Uo oc k^ + ^^ 4+ ^ and to different submm wavelengths 
using S Vo oc fo + ^ (see Equation IA9|) . 



